%This code produces figures 1,3,5

clear;
n=1;

rho=0.05;
mu_H=0.4;
mu_L=0.3;

M=1;
GG=mu_H/(rho+mu_H)*mu_L/(rho+mu_L);

up_bound_phi_L=-(mu_L+mu_H)/rho*M;

cutoff=1+2;% Change 2 to reflect the cost of external finance ($\kappa=2$ corresponds to the right-most graph of Figure 5) 
for phi_L=M:-0.0001:up_bound_phi_L+0.0001
if phi_L<0
phi_H=(M-mu_L*phi_L/(mu_L+mu_H))/mu_H*(mu_L+mu_H);

eta_1=(rho+mu_H)/phi_L;
eta_2=(rho+mu_L)/phi_H;


Phi_H=1/eta_2;
Phi_L=1/eta_1;

A=[eta_2 -mu_L/phi_H;-mu_H/phi_L eta_1];
om=eig(A);


om_1=om(2);
om_2=om(1);


c_1=1/(om_2-om_1)*om_2/om_1;
c_2=1/(om_1-om_2)*om_1/om_2;


x_m(n)=max(1/(om_2-om_1)*log(om_1^2/om_2^2*(eta_2-om_2)/(eta_2-om_1)),0);
x_m_norm(n)=x_m(n)./abs(phi_L);

der=(c_1*om_1*(1-om_1*Phi_H)*exp(-om_1*x_m(n))+c_2*om_2*(1-om_2*Phi_H)*exp(-om_2*x_m(n)))*(rho+mu_L)/mu_L;

derf=@(x) (c_1*om_1*(1-om_1*Phi_H)*exp(-om_1*x)+c_2*om_2*(1-om_2*Phi_H)*exp(-om_2*x))*(rho+mu_L)/mu_L -cutoff;


    if der>cutoff

        x_m_2(n)=fzero(derf,0.2*x_m(n));
        derf(x_m_2(n))
        
    else
        x_m_2(n)=x_m(n);

    end
else
    x_m(n)=0;
end
phi(n)=phi_L;
n=n+1;
end



xx=[M*sqrt(mu_L/mu_H) sqrt(mu_L)/sqrt(mu_H)*abs(phi-M)]./M;
zz1=[0 x_m];
zz2=[0 x_m_2];
zz=zz1-zz2;

yy1=[0 x_m];
yy2=[0 x_m_2];

FSsize=14;
sizedot=0.2;
ttarrow=0.1;



% Figure 1 in the text
figure
plot(xx,[0 x_m],'Linewidth',3)
box('off')
xlabel('coefficient of variation, $\theta$','Interpreter','latex','FontSize',FSsize)
ylabel('$X^\ast $','Interpreter','latex','FontSize',FSsize)



old_y_lim=ylim;
old_x_lim=xlim;

text(xx(min(find(zz1>0)))-2*sizedot, sizedot, 'A', 'interpreter','Latex','FontSize',FSsize)
text(xx(max(find(zz1>0)))-2*sizedot, sizedot, 'B', 'interpreter','Latex','FontSize',FSsize)
text(xx(end)-2*sizedot, sizedot, 'C', 'interpreter','Latex','FontSize',FSsize)


a=rectangle('Position',[xx(min(find(zz1>0)))-sizedot/2 -sizedot/2 sizedot sizedot],'FaceColor','red','Curvature',[1 1])
b=rectangle('Position',[xx(max(find(zz1>0)))-sizedot/2 -sizedot/2 sizedot sizedot],'FaceColor','red','Curvature',[1 1])
c=rectangle('Position',[xx(end)-sizedot/2 -sizedot/2 sizedot sizedot],'FaceColor','red','Curvature',[1 1])
ylim([0 old_y_lim(2)])
set (a, 'Clipping', 'off');
set (b, 'Clipping', 'off');
set (c, 'Clipping', 'off');

annotation('textarrow',[xx(min(find(zz1>0)))-ttarrow xx(min(find(zz1>0)))+ttarrow]/old_x_lim(2),[max(yy1)/2+ttarrow max(yy1)/2]/old_y_lim(2), 'String','$\theta_A \equiv \sqrt{\frac{\mu^L}{\mu^H}}$','Interpreter','latex','FontSize',FSsize)
annotation('textarrow',[xx(max(find(zz1>0)))+ttarrow xx(max(find(zz1>0)))-ttarrow]/old_x_lim(2),[max(yy2)/2.5 max(yy2)/2+ttarrow]/old_y_lim(2), 'String','$\theta_B \equiv (1-\lambda) \theta_A + \lambda \theta_C$','Interpreter','latex','FontSize',FSsize)
annotation('textarrow',[xx(min(find(zz1>0)))-ttarrow xx(min(find(zz1>0)))+ttarrow]/old_x_lim(2),[max(yy1)/3+ttarrow max(yy1)/2]/old_y_lim(2), 'String','$\theta_C \equiv \left(1+\frac{\mu^L+\mu^H}{\rho}\right) \sqrt{\frac{\mu^L}{\mu^H}}$','Interpreter','latex','FontSize',FSsize)

annotation('textbox',[0.4 0.7 0.6 0.08], 'String','Point A: $\phi^L=0                                                          $','Interpreter','latex','FontSize',FSsize,'LineStyle','none','FitBoxToText','off')
annotation('textbox',[0.4 0.6 0.6 0.08], 'String','Point C: $E_t\{ \int_t^\infty \phi_s e^{-\rho(t-s)} ds | \phi_t=\phi^L \} =0$','Interpreter','latex','FontSize',FSsize,'LineStyle','none','FitBoxToText','off')


% Figure 3

x_m_norm(1:min(find(zz1>0)))=NaN;
sizedot=0.2;
figure
plot(xx,[0 x_m_norm],'Linewidth',3)
box('off')
xlabel('coefficient of variation, $\theta$','Interpreter','latex','FontSize',FSsize)
ylabel('$\frac{X^\ast}{|\phi_L|} $','Interpreter','latex','FontSize',FSsize,'Rotation',0)

xline(xx(min(find(zz1>0)))-0.5*sizedot,'--')


old_y_lim=ylim;
old_x_lim=xlim;
text(xx(min(find(zz1>0)))-2*sizedot, sizedot, '$\theta_A$', 'interpreter','Latex','FontSize',FSsize)
text(xx(max(find(zz1>0)))-2*sizedot, sizedot, '$\theta_B$', 'interpreter','Latex','FontSize',FSsize)
text(xx(end)-2*sizedot, sizedot, '$\theta_C$', 'interpreter','Latex','FontSize',FSsize)


a1=rectangle('Position',[xx(min(find(zz1>0)))-sizedot -1.5*sizedot sizedot 3*sizedot],'FaceColor','red','Curvature',[1 1])
b1=rectangle('Position',[xx(max(find(zz1>0)))-sizedot/2 -1.5*sizedot sizedot 3*sizedot],'FaceColor','red','Curvature',[1 1])
c1=rectangle('Position',[xx(end)-sizedot/2 -1.5*sizedot sizedot 3*sizedot],'FaceColor','red','Curvature',[1 1])
ylim([0 old_y_lim(2)])
set (a1, 'Clipping', 'off');
set (b1, 'Clipping', 'off');
set (c1, 'Clipping', 'off');



% Figure 5
sizedot=0.2;
figure
plot(xx,[0 x_m],'-r',xx,[0 x_m_2],'--b','Linewidth',3)
box('off')
xlabel('coefficient of variation of $\theta$','Interpreter','latex','FontSize',FSsize)
ylabel('$X^\ast, X^{\ast \ast}$','Interpreter','latex','FontSize',FSsize,'Rotation',0)

xline(xx(max(find(zz>0))),'-.')
old_y_lim=ylim;
old_x_lim=xlim;

text(xx(min(find(zz1>0))), sizedot, '$\theta_A$', 'interpreter','Latex','FontSize',FSsize)
text(xx(max(find(zz1>0))), sizedot, '$\theta_B$', 'interpreter','Latex','FontSize',FSsize)
text(xx(end), sizedot, '$\theta_C$', 'interpreter','Latex','FontSize',FSsize)
text(xx(max(find(zz>0))),sizedot,'$\theta_D$', 'interpreter','Latex','FontSize',FSsize)


a=rectangle('Position',[xx(min(find(zz1>0)))-sizedot/2 -sizedot/2 sizedot sizedot],'FaceColor','red','Curvature',[1 1])
b=rectangle('Position',[xx(max(find(zz1>0)))-sizedot/2 -sizedot/2 sizedot sizedot],'FaceColor','red','Curvature',[1 1])
c=rectangle('Position',[xx(end)-sizedot/2 -sizedot/2 sizedot sizedot],'FaceColor','red','Curvature',[1 1])
d=rectangle('Position',[xx(max(find(zz>0)))-sizedot/2 -sizedot/2 sizedot sizedot],'FaceColor','red','Curvature',[1 1])
ylim([0 old_y_lim(2)])
set (a, 'Clipping', 'off');
set (b, 'Clipping', 'off');
set (c, 'Clipping', 'off');
set (d, 'Clipping', 'off');

annotation('textarrow',[xx(min(find(zz1>0)))-ttarrow xx(min(find(zz1>0)))+ttarrow]/old_x_lim(2),[max(yy1)/2+ttarrow max(yy1)/2]/old_y_lim(2), 'String','$X^{\ast}$','Interpreter','latex','FontSize',FSsize)
annotation('textarrow',[xx(max(find(zz1>0)))+ttarrow xx(max(find(zz1>0)))-ttarrow]/old_x_lim(2),[max(yy2)/2 max(yy2)/2+ttarrow]/old_y_lim(2), 'String','$X^{\ast \ast}$','Interpreter','latex','FontSize',FSsize)

